clear all
close all

CaV = 1; %Ca: -1, Voltage:1
legendtext = 'ms';
clim = [0 250]; % color lim 0-250msec

pointX = 128; %reference X
pointY = 128; %reference Y

tStart = 1;
tStop = 1999; %usually one less than number of frames

MinPeakProminence = 0.007; %parameter define for picking peaks at reference 
                         %(depending on signal to noise ratio): 
                         %default  = 0.01 1/14/2019 GSD data original
                         %default  = 0.007 1/14/2019 GSD data diff
                         
peakWindow = [-150 250]; %ms window for picking peaks correlated with reference peak

sizeFOV = 16.5647; %1x lens:16.5647 and 1.6x lens:10.35 mm
%80mm is the backfocal distance for the Leica Planapo lenses (0.63x, 1.0x, 1.6x, 2.0x).  
%The focal length of the projection lens / the backfocal distance of the objective lens * 
%the objective lens magnification = system magnification.  
%In the case of the 1.6x lens as objective with 85mm lens as projector:  
%85/80*1.6 = 1.7x.  The N256 sensor size is 17.6x17.6, so 17.6/1.7=10.35mm field of view
%In the case of the 1.0x lens as objective with 85mm lens as projector:  
%85/80*1.0 = 1.0625x.  The N256 sensor size is 17.6x17.6, so 17.6/1.0625=16.5647mm field of view

% field of view of interest for activation speed calculation
windowSpeedX = [61 200]; %pixel
windowSpeedY = [61 200]; %pixel

% size of window for local activation speed calculation
windowLocalSpeedX = 11; %pixel
windowLocalSpeedY = 11; %pixel

actSpeedMax = 50; %cm/sec
actSpeedMin = 0; %cm/sec



%get filename and pathname of GSD file
[FileName,PathName,FileIndex] = uigetfile('*.gsd','Select the binary OMS file', ... 
    'C:\Users\spark71\Documents\MATLAB\OMS\V2025_test');

filenameInput = [PathName FileName];





%read read oms GSD file into mat file
out = readoriginoms2020(filenameInput, CaV);

xsize = out.xsize; %size of total pixel at x-axis 
ysize = out.ysize; %size of total pixel at y-axis 
frameperSec = out.fs; %Frames/sec
umperPixel = sizeFOV/xsize*1000; % um
totalframenumber = out.totalframenumber; %total frame number

rawData = out.data; %[x y t]
bgImageData = out.bgimagedata; 
diffImageData = out.differentialdata;


%filtering data using time and space filters
Data_t_filtered = smoothers2020(rawData, 'time', 4);
Data_ts_filtered = smoothers2020(Data_t_filtered, 'space', 5, 4);

Data_anal = Data_ts_filtered; 
rawData_anal = rawData;

% Data_anal = diff(Data_ts_filtered,1,3); 
% rawData_anal = diff(rawData,1,3);






% find peaks/peak locations at reference (pointX, pointY)

dataPt = double(squeeze(Data_anal(pointX,pointY,tStart:tStop)));
rawDataPt = double(squeeze(rawData_anal(pointX,pointY,tStart:tStop)));

[pksPt, locsPt] = findpeaks(dataPt, 1:length(dataPt), 'MinPeakProminence', MinPeakProminence); % reference peak location

hff = figure;
plot((1:length(dataPt))*(1000/frameperSec), rawDataPt)
hold on
plot((1:length(dataPt))*(1000/frameperSec), dataPt)
plot(locsPt*(1000/frameperSec), pksPt, '+')
xlabel('time(ms)');
ylabel('signal');
box off
%saveas(hff, [outFileName '_signal'], 'fig')


% find peaks/peak locations over all pixels

pks =  zeros(xsize, ysize, 100);
locs =  zeros(xsize, ysize, 100);
pksN = zeros(xsize, ysize); % pk number at i and j
for i = 1:xsize
    i
    for j=1:ysize
        
        dataIJ = double(squeeze(Data_anal(i,j,tStart:tStop)));
        [pksIJ, locsIJ] = findpeaks(dataIJ, 1:length(dataIJ), 'MinPeakProminence', MinPeakProminence);
    
        if(isnan(pksIJ))
        else
            pksN(i,j) = length(pksIJ);
            if(length(pksIJ)<size(pks,3))
                for k = 1:length(pksIJ)
                    pks(i,j,k) = pksIJ(k);
                    locs(i,j,k) = locsIJ(k);
                end
            else
                for k = 1:size(pks,3)
                    pks(i,j,k) = pksIJ(k);
                    locs(i,j,k) = locsIJ(k);
                end
            end
        end 
    end
end




% pick only activation time coordinated with reference peaks and remove all
% other peaks using function mappingImage2020A with maxlevel and minlevel
numberOfPulses = length(locsPt); % number of peak at reference;
for pp = 1:numberOfPulses

    inputData = (locs(:,:,pp)')*(1000/frameperSec); %msec

    outputfilename = [FileName(1:(length(FileName)-4)) '_p' num2str(pp)];
    
    %clim = [-50 200];
    maxlevel = locsPt(pp)*(1000/frameperSec)+peakWindow(2);
    minlevel = locsPt(pp)*(1000/frameperSec)+peakWindow(1);

    outMapping(pp) = mappingImage2020A(inputData, xsize, ysize, outputfilename, legendtext, clim, maxlevel, minlevel);
end

% average all peak's activation time
actTimeTemp = zeros(xsize, ysize);
for pp = 1:numberOfPulses
    actTimeTemp = actTimeTemp + outMapping(pp).zim; %msec
     
    if(pp == numberOfPulses)
        actTimeAvg = actTimeTemp/numberOfPulses;
        reAct = reshape (actTimeAvg, [xsize*ysize 1]);
        tempIndex = find(reAct>=0);
        maxTime = max(reAct(tempIndex));
        minTime = min(reAct(tempIndex));
        actTimeAvgAdjusted = actTimeAvg - minTime;
    end
end

maxlevel = maxTime+10;
minlevel = minTime-10;
    
outAvg = mappingImage2020A(actTimeAvgAdjusted, xsize, ysize, [FileName(1:(length(FileName)-4)) '_Avg'], legendtext, clim, maxlevel, minlevel);



%calculate the global activation speed

meanActTimeX = mean(actTimeAvgAdjusted(windowSpeedX(1):windowSpeedX(2), windowSpeedY(1):windowSpeedY(2)),1);
meanActTimeY = mean(actTimeAvgAdjusted(windowSpeedX(1):windowSpeedX(2), windowSpeedY(1):windowSpeedY(2)),2)';

meanActTimeX0 = min(meanActTimeX);
meanActTimeY0 = min(meanActTimeY);

meanActTimeXX = meanActTimeX - meanActTimeX0; %msec
meanActTimeYY = meanActTimeY - meanActTimeY0; %msec

travelLocationX = umperPixel*(1:length(meanActTimeXX)); %um
travelLocationY = umperPixel*(1:length(meanActTimeYY)); %um

linearFittingX = polyfit(travelLocationX, meanActTimeXX, 1); %msec/um ->sec/mm
linearFittingY = polyfit(travelLocationY, meanActTimeYY, 1); %msec/um ->sec/mm

actSpeedGlobal = 0.1/(linearFittingX(1)^2+linearFittingY(1)^2)^0.5; % cm/sec
actSpeedX = 10*linearFittingX(1)*actSpeedGlobal; % NonDimensional Unit
actSpeedY = 10*linearFittingY(1)*actSpeedGlobal; % NonDimensional Unit
actSpeedAngle = atan2(actSpeedX, actSpeedY); %rad
actSpeedAngle2 = atan2(actSpeedX, actSpeedY)/pi*180; %degree
actSpeedX2 = actSpeedGlobal*actSpeedX; % cm/sec
actSpeedY2 = actSpeedGlobal*actSpeedY; % cm/sec


hff2 = figure;
plot(travelLocationX, meanActTimeXX)
hold on
plot(travelLocationX, linearFittingX(1)*travelLocationX+linearFittingX(2))
xlabel('location(um)');
ylabel('time(ms)');
box off

hff3 = figure;
plot(travelLocationY, meanActTimeYY)
hold on
plot(travelLocationY, linearFittingY(1)*travelLocationY+linearFittingY(2))
xlabel('location(um)');
ylabel('time(ms)');
box off



%calculate the local activation speed

actSpeedLocal = zeros((windowSpeedX(2)-windowSpeedX(1)+1),(windowSpeedY(2)-windowSpeedY(1)+1)); % cm/sec
actSpeedLocalX = zeros((windowSpeedX(2)-windowSpeedX(1)+1),(windowSpeedY(2)-windowSpeedY(1)+1)); % NonDimensional Unit
actSpeedLocalY = zeros((windowSpeedX(2)-windowSpeedX(1)+1),(windowSpeedY(2)-windowSpeedY(1)+1)); % NonDimensional Unit
actSpeedLocalAngle = zeros((windowSpeedX(2)-windowSpeedX(1)+1),(windowSpeedY(2)-windowSpeedY(1)+1)); %rad
actSpeedLocalAngle2 = zeros((windowSpeedX(2)-windowSpeedX(1)+1),(windowSpeedY(2)-windowSpeedY(1)+1)); %degree
actSpeedLocalX2 = zeros((windowSpeedX(2)-windowSpeedX(1)+1),(windowSpeedY(2)-windowSpeedY(1)+1)); % cm/sec
actSpeedLocalY2 = zeros((windowSpeedX(2)-windowSpeedX(1)+1),(windowSpeedY(2)-windowSpeedY(1)+1)); % cm/sec

for pp = 1:(windowSpeedX(2)-windowSpeedX(1)+1)
    pp
    xIndex = windowSpeedX(1)+pp-1;
    for qq = 1:(windowSpeedY(2)-windowSpeedY(1)+1)
        yIndex = windowSpeedY(1)+qq-1;
        
        actTimeLocalInterest = actTimeAvgAdjusted((xIndex-round((windowLocalSpeedX-1)/2)):(xIndex+round((windowLocalSpeedX-1)/2)), (yIndex-round((windowLocalSpeedY-1)/2)):(yIndex+round((windowLocalSpeedY-1)/2)));
        meanActTimeLocalX = mean(actTimeLocalInterest,1);
        meanActTimeLocalY = mean(actTimeLocalInterest,2)';
        
        meanActTimeLocalX0 = min(meanActTimeLocalX);
        meanActTimeLocalY0 = min(meanActTimeLocalY);

        meanActTimeLocalXX = meanActTimeLocalX - meanActTimeLocalX0; %msec
        meanActTimeLocalYY = meanActTimeLocalY - meanActTimeLocalY0; %msec

        travelLocationLocalX = umperPixel*(1:length(meanActTimeLocalXX)); %um
        travelLocationLocalY = umperPixel*(1:length(meanActTimeLocalYY)); %um

        linearFittingLocalX = polyfit(travelLocationLocalX, meanActTimeLocalXX, 1); %msec/um ->sec/mm
        linearFittingLocalY = polyfit(travelLocationLocalY, meanActTimeLocalYY, 1); %msec/um ->sec/mm

        actSpeedLocal(pp,qq) = 0.1/(linearFittingLocalX(1)^2+linearFittingLocalY(1)^2)^0.5; % cm/sec
        actSpeedLocalX(pp,qq) = 10*linearFittingLocalX(1)*actSpeedLocal(pp,qq); % NonDimensional Unit
        actSpeedLocalY(pp,qq) = 10*linearFittingLocalY(1)*actSpeedLocal(pp,qq); % NonDimensional Unit
        actSpeedLocalAngle(pp,qq) = atan2(actSpeedLocalX(pp,qq), actSpeedLocalY(pp,qq)); %rad
        actSpeedLocalAngle2(pp,qq) = atan2(actSpeedLocalX(pp,qq), actSpeedLocalY(pp,qq))/pi*180; %degree
        actSpeedLocalX2(pp,qq) = actSpeedLocal(pp,qq)*actSpeedLocalX(pp,qq); % cm/sec
        actSpeedLocalY2(pp,qq) = actSpeedLocal(pp,qq)*actSpeedLocalY(pp,qq); % cm/sec
        
    end
end

dataHistoActSpeedLocal = reshape(actSpeedLocal,1,(windowSpeedX(2)-windowSpeedX(1)+1)*(windowSpeedY(2)-windowSpeedY(1)+1));
dataHistoActSpeedLocalAngle = reshape(actSpeedLocalAngle,1,(windowSpeedX(2)-windowSpeedX(1)+1)*(windowSpeedY(2)-windowSpeedY(1)+1));

isINF1 = isinf(dataHistoActSpeedLocal);
dataINF = find(isINF1==1);
dataHistoActSpeedLocal(dataINF) = NaN;
dataOutRange1 = find(dataHistoActSpeedLocal>actSpeedMax);
dataHistoActSpeedLocal(dataOutRange1) = NaN;
dataOutRange2 = find(dataHistoActSpeedLocal<actSpeedMin);
dataHistoActSpeedLocal(dataOutRange2) = NaN;

isNAN1 = isnan(dataHistoActSpeedLocal);
actSpeedLocalDataN = length(dataHistoActSpeedLocal)-sum(isNAN1);
[x, y, p, q] =  normfit(dataHistoActSpeedLocal(find(isNAN1==0)));
actSpeedLocalMean = x;
actSpeedLocalStd2 = y;
actSpeedLocalStdErr = y/(actSpeedLocalDataN^0.5);
actSpeedLocalMeanLimit = p;
actSpeedLocalStdLimit = q;
actSpeedLocalStdErrLimit = q/(actSpeedLocalDataN^0.5);

isINF2 = isinf(dataHistoActSpeedLocalAngle);
actSpeedLocalAngleDataN = length(dataHistoActSpeedLocalAngle)-sum(isINF2);
[x2, y2, p2, q2] =  normfit(dataHistoActSpeedLocalAngle(find(isINF2==0)));
actSpeedLocalAngleMean = x2;
actSpeedLocalAngleStd2 = y2;
actSpeedLocalAngleStdErr = y2/(actSpeedLocalAngleDataN^0.5);
actSpeedLocalAngleMeanLimit = p2;
actSpeedLocalAngleStdLimit = q2;
actSpeedLocalAngleStdErrLimit = q2/(actSpeedLocalAngleDataN^0.5);

hff4 = figure;
hist(dataHistoActSpeedLocal(find(isNAN1==0)),actSpeedMin:1:actSpeedMax)
box off
xlabel('Speed (cm/sec)', 'fontsize', 14)
ylabel('Count', 'fontsize', 14)
xlim([actSpeedMin actSpeedMax])

hff5 = figure;
hist(dataHistoActSpeedLocalAngle(find(isINF2==0))/pi*180,-180:5:180)
box off
xlabel('Angle(degree)', 'fontsize', 14)
ylabel('Count', 'fontsize', 14)
xlim([-180 180])



%calculate the activation speed of tissue at each pulse

actTimePulse = zeros(xsize, ysize, numberOfPulses);
for pp = 1:numberOfPulses

    actTimePulseTemp = outMapping(pp).zim; %msec
    reActTimePulseTemp = reshape (actTimePulseTemp, [xsize*ysize 1]);
    tempPulseIndex = find(reActTimePulseTemp>=0);
    minPulseTime = min(reActTimePulseTemp(tempPulseIndex));
    actTimePulseTempAdjusted = actTimePulseTemp - minPulseTime;
    actTimePulse(:,:,pp) = actTimePulseTempAdjusted; %msec

end
